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ACCURACY STUDY OF FINITE DIFFERENCE METHODS 


By Nancy Jane Cyrus and Robert E. Fulton 
Langley Research Center 

SUMMARY 

A method for studying the accuracy of finite difference approximations for linear 
differential equations is presented and utilized. Definitive expressions for the error in 
each approximation are obtained by using Taylor series to derive the differential equa- 
tions which exactly represent the finite difference approximations. The resulting differ- 
ential equations are accurately solved by a perturbation technique which yields the error 
directly. 

This method is used to assess the accuracy of two alternate forms of central finite 
difference approximations for solving boundary value problems in structural analysis 
which are governed by certain equations containing variable coefficients. A "half station 
approximation" in which finite difference approximations are made before expanding 
derivatives of function products is compared with a "whole station approximation" in which, 
derivatives of function products are expanded first for string, beam, and axisymmetric 
circular plate problems. An example of a square membrane is given as an application of 
the method to partial differential equations. 

INTRODUCTION 

The differential equations governing the behavior of structural boundary value prob- 
lems are often solved by approximating the derivatives by finite differences and solving 
the resulting algebraic equations on a digital computer. For complicated structures the 
number of simultaneous equations resulting from finite difference approximations can be 
sufficiently large to exceed the capacity of the computer or introduce round-off error. 

For such problems, the accuracy of the difference procedure can be a critical item in 
obtaining meaningful design results. In reference 1, for example, it was found that accu- 
rate answers for the stress in a shell could not be obtained by using certain finite differ- 
ence approximations unless the mesh spacing was smaller than machine capacity 
permitted. 

The most popular difference approximations are the so-called central differences 
which are given in textbooks on numerical methods. There are alternate formulations of 
central differences which can be used when odd order derivatives occur in the differential 


equations. Such a situation exists in structural problems, for example, where inplane 
loads are not uniform (a column loaded by its own weight or a shell of revolution sub- 
jected to arbitrary loads) or where the stiffness of the structure is nonuniform (a tapered 
beam or a variable thickness shell). 

In this paper a method for studying the accuracy of finite difference approximations 
is presented and utilized. As illustrative examples, the method is used to assess the 
accuracy of two alternate forms of central finite difference approximations used in struc- 
tural problems through application to string, beam, axisymmetric circular plate, and 
square membrane problems. The same approach can be used to evaluatie the accuracy of 
finite element methods. 


SYMBOLS 

B(y),L(y) linear differential operators 
D(y),E(y) linear difference operators 
f(x) nondimens ional tension in a beam or string 

g(x) nondimensional stiffness of beam 

h finite difference spacing 

i any integer 

k superscript describing set of boundary conditions 

m,n Fourier wave numbers 

p(x) nondimensional lateral load 

q(x) boundary condition 

x,z descriptive coordinates of beam, string, plate, or membrane 

y deflection of beam, string, plate, or membrane 

Y deflection function in perturbation series 


2 


slope of plate 


' 4 > 

r boundary curve 

Prime or Roman numeral with a symbol denotes differentiation with respect to x. 

METHOD OF APPROACH 

The usual approach in a finite difference accuracy study is to carry out the numeri- 
cal solution to a number of problems for which exact solutions can be obtained and to com- 
pare the resulting numerical answers at each station with the exact answers. Such a 
procedure has the liability that comparisons can only be made for each problem at spe- 
cific stations and the calculations must be redone each time the mesh size changes. 

The approach used in this paper is to isolate the principal finite difference error so 
that its magnitude and character can be evaluated. The finite difference approximations 
are expanded in Taylor series to give differential equations which are exactly equivalent 
to the finite difference approximations. Solving the resulting differential equations by a 
perturbation technique yields analytical expressions for the principal error term. These 
expressions are independent of mesh spacing and give a clear indication of the accuracy 
of the difference approximations not just at discrete points but over the whole domain of 
interest. 

Consider the differential equation 


L(y) = p (1) 

with a necessary and sufficient set of k boundary conditions, each of the form 

B k (y) = q k (on T) (2) 


Equation (1) may represent either an ordinary or partial differential equation. For 
example, equation (1) takes the form for a string of 

dx 2 

and for a membrane of 

V 2 y = p(x,z) 
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where V 2 is the Laplacian operator 


V 2 = -i£- + 

9x 2 8z 2 

The differential equation (1) is approximated by finite differences and is replaced by 
a finite difference recursion formula at the ith station of the form 


D (Yi) = Pi 


( 3 ) 


where D(yi) is the equivalent finite difference operator for L(y) and is expressed in 
terms of y evaluated at the appropriate finite difference stations. 

A similar treatment for each of the k boundary conditions leads to replacing 
equation (2) by 


E k ( yi ) = qi k (on T) (4) 


where E k is the finite difference operator for B k . Note that operators of the form of 
equations (3) and (4) also result for finite element problems if a continuum is approxi- 
mated by finite elements and the approximate equilibrium equations and boundary condi- 
tions are obtained. 


The finite difference recursion equations (3) and (4) may be expanded about the ith 
point by using the appropriate Taylor series expansion, such as the one -dimensional 
expansion 


y i± l = Yi ± h Yi ’ + 


h 2 yj > 

2 '. 


h 3 


yi 


hV v 

4! 


For any central finite difference method, the order of error of the approximation is 
proportional to h 2 and the finite difference recursion formula takes the form 


D (yi) = L o(yi) + h2L i(yi) + h4L 2(yi) + • • • 


(5a) 


where Lq, Lj, and L 2 are differential operators which depend on the approximation 
method used. A similar treatment for E k leads to 


E k ( Yi ) = B 0 k (yi) + h 2 B 1 k (y i ) + h 4 B 2 k (yi) + . . . (5b) 


For other difference approximations all powers of h may occur. 
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By use of equations (5a) and (5b), the finite difference equations (3) and (4) now 
take the form 

L o(yi) - Pi + h2l n(yi) + h4L 2(yi) + • • • - 0 (6) 

B o k (yi) - q i k + h2B i k (yi) + h4 B 2 k (yi) + • • • = 0 (on d (7) 

Equation (6) and its k boundary condition equation (7) are the differential equations 
which represent the finite difference recursion equations (3) and (4). As the increment 
h goes to zero, equation (6) and equation (7) should approach equation (1) and equa- 
tion (2), respectively. In fact, if the finite difference approximations used are a con- 
vergent set, then 

L 0 = L 

and 

B 0 = B 

The solution to equations (6) and (7) gives an analytical representation of the numer- 
ical finite difference answers. Unfortunately, because of the infinite number of terms in 
equation (6) a closed form solution does not appear feasible. However, in a practical 
problem where the size of the region is scaled to be of the order one, h is perhaps 0.1 
or 0.01 or even smaller. This small value of h suggests that equation (6) may be 
solved by a perturbation technique with the perturbation parameter taken to be h^. 

Let the solution to equations (6) and (7) be taken in the form 

Yi = Y 0 + h 2 Y x + h 4 Y 2 + . . . (8) 

Substituting equation (8) into equation (6) leads to 

L o( Y o) - Pi + h2 [ L o( Y l) + L l( Y o)] + h4 [ L o( Y 2 ) + L l( Y l) + L 2 (Y 0 )] + . . . = o (9) 

subject to k boundary conditions of the form 

B oVo) - <li k + h2 [ B o k ( Y l) + B l k ( Y o)] + h 4 [ B 0 k ( Y 2) + B l k ( Y l) + B 2 k ( Y o)] +•••=» 

( 10 ) 

If each order of error term is solved in sequence, the following series of problems 
result: 

V Y o) - Pi ■ 0 ( B O k ( Y o) - «i k - °) ( 11 ) 

M Y l) + L l( Y o) = 0 ( B 0 k ( Y l) + B l k Y o) * o) (12) 

l o( Y 2) + L l( Y l) + L a( Y o) - 0 ( B 0 k ( Y 2) + B l k ( Y l) + B 2 k ( Y o) - °) ( 13 ) 

(...) 
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If the finite difference approximation is a convergent one, equation (11) is equation (1) and 
Yp given by equation (11) is the exact solution to equation (1). 

From the form of in equation (8) it is seen that Yj can be interpreted as the 
principal error term in the finite difference results in relation to the exact answer to the 
problem. If two different finite difference approximations are to be considered, a com- 
parison of error terms Yj resulting from the two different approximations indicates the 
relative accuracy of the two approximations. 

ILLUSTRATIVE PROBLEMS 


To illustrate the procedure, the results from two finite difference methods are com- 
pared for a string with nonuniform tension. Several additional structural problems are 
given in the appendices as further examples of the use of the method. The examples 
include a beam with nonuniform stiffness, an axisymmetric plate, and a square membrane 
subjected to a sinusoidal loading. The beam and string problems were taken from refer- 
ence 2 and are given for completeness. 


String With Nonuniform Tension 

Consider a string of constant length with nonuniform tension f(x) subjected to a 
lateral load p(x). The governing differential equation is 

(fy')' + p(x) = 0 (14) 

where primes denote differentiation with respect to x. The variables are nondimension- 
alized so that the length of the string is one and the tension at the left end is one. The 
boundary conditions are 

y(xp) = 0 y(x 0 + l) = 0 (15) 

where xq is the coordinate at the left end of the string. Equation (14) may be solved 
by dividing the string into stations of equal spacing h. The finite difference equations 
are written in terms of displacements at the ith station (i = 1, 2, 3, . . .). 

Two finite difference approximations to be considered are denoted, for convenience, 
the "half station" approximation and the "whole station” approximation. For equa- 
tion (14) these two approximations lead to the following finite difference expressions: 


Half station approximation 


(fy')i’ +P i 




+ p i 




+ Pi =0 


(16) 
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II 


or whole station approximation 


f. f. * 

(fy" + f’y')i + Pi = 4(yn - 2y A + y i+1 ) + + y i+1 ) + Pj 




f i - ~§“Vi-i - 2f i y i + ( f i + ^-k+i 


+ Pi = 0 


( 17 ) 


Note that the half station approximation is the natural result of making the finite differ- 
ence approximation before expanding the derivatives, whereas the whole station approxi- 
mation results from making the approximation after the expansion. The latter type of 
approximation is widely used. (See, for example, refs. 3 and 4.) Both of the preceding 

p 

sets of finite difference approximations can be shown to be of order h and yet they 
clearly lead to different coefficients for the simultaneous equations in terms of the dis- 
placements at the ith station. Of concern here are the relative magnitudes of the errors 
in these two different approximations. 

Expand the finite difference recursion equations (16) and (17) about the ith point by 
using such Taylor series expansions as 

h 2 

y i±l =y i * hy i' + 2T y i" * ’ ‘ • 


f i ± l 


fi ± hfj 


*2 

+ S-f i 

21 1 


For both the half station and whole station approximations, this procedure leads to a dif- 
ferential equation of the form given by equation (6) where 


L 0( y i)= ( f i y i')’ 
B o(yi) = y i 


(18) 


and for the half station approximation 

fm IV W" W' 


L i( yi ) 12 

B l( y i) = 0 


8 24 




(19) 


fi7i VI . fi’yi v . f i" y i IV . fryf + VV + fiV 


L2 ( y i) 360 ' 120 

B 2(yi) = o 


96 


144 384 1920 




( 20 ) 
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and for the whole station approximation 


L l( y i) 

B l( y i) 

L 2(yi) 

B 2(yi) 


W v W" 
+ 


12 


6 


kvi" fi’yi 1 


360 


120 


( 21 ) 


( 22 ) 


Equations (6) and (7) with equations (18) and either equations (19), (20), . . ., or 
equations (21), (22), . . ., are clearly differential equations and associated boundary con- 
ditions which represent exactly the two finite difference recursion equations (16) and (17) 
and their associated boundary conditions. 

By using the method described in the previous section, the principal error func- 
tions Yj defined by equation (8) corresponding to the half station and the whole station 
finite difference approximations have been obtained for a family of problems. These 
problems are a string having a lateral load which is distributed uniformly and a tension 


force f(x) which varies as follows: 

f(x) = i 
x 1A 

(1 £ n £ 6) 

subject to the boundary conditions 

y(i) = o 



y(2) = o 


and 

f (x) = 1 + X 11 

(2 2 n ^ 6) 

subject to the boundary conditions 

y(0) = 0 



y(i) = o 



Where f(x) is linear (corresponding to f(x) = 1, x, or 1 + x), the results for the 
half station and whole station finite difference approximations are exactly the same. In 
fact, for f(x) = 1, both difference answers are the exact answer. For all other cases, 
however, the two difference methods lead to different results. It is useful to compare 
the results for f(x) = in detail as a typical example. 

X'* 
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( 23 ) 


£ . 31 
5 75 


For f(x) = \ and y(l) = y(2) = 0, 

x 6 

y 0 = -; 

and the half station approximation is 

Yl = - 41 x 4 + ^ 
Yl 1125 6 

and the whole station approximation is 

v = 187 -4 , 4 

Yl 450 3 


x 4 _ 16 


X 75 


31_ x 2 + 

86 

150 X + 

1125 

. 31x2. 

26 

30 

225 


(24) 


(25) 


The two error terms Yj over the length of the string are presented in fig- 
ure 1(a). Solutions were also obtained for the error terms in deflection for all the 
remaining load functions f(x) noted previously. Additional results for f(x) = 1 + x 8 
are shown in figure 1(b). The remaining solutions are not shown because figure 1 serves 
to illustrate the character of the results. An overall measure of the relative errors in 
the two methods is shown later for all solutions obtained. 


Although errors in the deflections of the string are important, errors in numerically 
obtained derivatives should also be considered for a thorough error analysis. Therefore, 
results were obtained by using the finite difference answers for approximate curvatures 
(second derivatives). The second difference operator was applied to the difference 
results; Taylor and perturbation series expansions were then applied to yield 


y i" - J2( y i-1 ' 2y i + y i + l) = Y °" + h V + 15 + ‘V 



or 


Y' = Y 0 " + h 2 (Y 1 " + ^) + . . . 


2 1 
The h error terms in the curvatures for f(x) = — ^ are as follows: 

x J 


For the half station approximation, 


IV 


Y t" + — — = - Ml x 2 - x + 

1 12 375 75 

and for the whole station approximation, 

v IV 

Y " + I0_-- 874 x 2 + 6x-M 
Y 1 + 12 " 75 X +bX 25 


(26) 


(27) 


(28) 
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The error in the curvature for each of the two approximations is also given in fig- 

1 n 

ure 1(a) for f(x) = — and in figure 1(b) for f(x) = 1 + x . Results for the remaining 
x3 

load functions are shown in a subsequent section in the form of an overall measure of the 
relative error. 

Numerical calculations were also carried out for the deflections and curvatures 
for the problems cited to determine whether the analytical errors adequately represented 
the numerical errors. The data are not included herein; however, for h less than about 
0.1 all analytical errors agree with calculated numerical errors within 1 percent. 

Beam, Plate, and Membrane Examples 

Appendix A contains examples of a simply supported beam having a nonuniform 
bending stiffness and subjected to a uniformly distributed load. Figure 2 shows the dis- 
tribution of deflection and curvature errors for a linearly tapered beam. Examples of a 
clamped circular plate and a simply supported annular plate under uniformly distributed 
load are given in appendix B. Appendix C contains results for a square membrane sub- 
jected to a single term Fourier load. 

RELATIVE ERRORS OF THE HALF STATION AND 
WHOLE STATION APPROXIMATIONS 

Although results such as those given in figures 1 and 2 are usually sufficient to 
identify which of the two approximations is superior for a given problem, identification 
of the superior method for specific results is sometimes difficult (for example, the 
curvature errors of fig. 1(b)). Moreover, a quantitative measure of the relative accuracy 
of the approximations is desirable. Probably the fairest comparison of their overall 
merit can be made by examining the root-mean- square values of the errors for the whole 
structure; that is, 

■ icr Y ‘ 2dx 

for the error in deflection and 



for the error in curvature, where the integration is over the (unit) length of the string, 
beam, or plate. Thus, to assess quantitatively the relative merits of the half station and 
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whole station approximations for the various problems solved, the ratios 

ov. , ,,Av and ov »» /(Tv " have been calculated for each problem. 

Y l,half/ Y 1>whole Yj half / Yj >whole 

The results are shown in figure 3. 

DISCUSSION OF RESULTS OF SAMPLE PROBLEMS 

The results given in figure 3(a) show that for all problems studied, the error in 
the deflection resulting from use of the half station approximation is less than the error 
obtained from use of the whole station approximation, in some problems by an order of 
magnitude. The investigation of the accuracy of the curvature approximations gives the 
same result in general. Thus, the half station method is usually superior for calculation 
of both deflection and bending curvature for the problems studied. 

Although the results clearly favor the half station approximation, one exception 
occurs: for the string with the load f(x) = 1 + x^, the error in the curvature is 25 per- 
cent greater with the half station approximation. The difference between the two approxi- 
mations is seen to be generally less in calculating the second derivatives of deflections 
than in calculating the deflections themselves. 

The analytical representation of errors in the present paper shows the danger of 
using numerical data at a single station or a few points to characterize the error in a 
problem. An example is shown in figure 1(a) for f(x) = -ij. If comparisons are made of 
the curvature near the end x = 1, the whole station approximation appears much more 
accurate than the half station approximation; whereas figure 3(b) shows clearly that the 
average error with the whole station approximation is over twice as great. 

The present approach to error assessment may also be useful for comparison of 
different finite element structural approximations. In fact, the recursion formulas 
given by the half station approximation (eq. (16) and eq. (A2)) are the same recursion 
formulas that occur for a finite element model consisting of rigid bars connected by 
rotational springs, which often is used to represent a physical problem such as a beam- 
column (for example, ref. 5). Thus, the results of the present study verify that the finite 
element model of reference 5 is a good representation of the behavior of the continuum 
problem. 

A practical consideration which supports the use of the half station method is the 
symmetry of the matrix of coefficients in this approximation. By contrast, the matrix 
of coefficients associated with whole stations is not symmetric. Matrix symmetry can 
be of great value for many numerical procedures associated with eigenvalue routines and 
simultaneous equation solving routines and, in some problems, is required for an effi- 
cient numerical solution of a large order system. 
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The results for the square membrane example given in appendix C demonstrate the 
application of the method to partial differential equations and indicate the relative accu- 
racy of two alternate patterns for the Laplacian operator. The conventional pattern 
having error of order h 2 is compared with a so-called refined pattern which can be 
shown to have an error of order h 4 if the Laplacian of the loading vanishes. It is seen 
that for a single Fourier loading the standard pattern is actually more accurate than the 
refined pattern. Definitive expressions for the error terms are presented for both 
approximations. These expressions give the number of finite difference stations which 
are required within the length of a deflection Fourier wave to restrict finite difference 
answers to a given percentage of error. 

CONCLUDING REMARKS 

A new procedure has been developed to determine an analytical representation of 
the error in a finite difference solution to a specified problem. This procedure allows 
a direct comparison, independent of mesh size, between difference approximations. The 
procedure appears to have considerable merit for assessment of the relative accuracy of 
finite difference and finite element numerical techniques of linear structural analysis. 

By using this procedure, a comparison has been made of the accuracy of two finite 
difference approximations for solving structural problems through applications to a 
spectrum of beam and string problems having the characteristics of nonuniform stiff- 
ness and inplane load and to two circular plate problems. The methods investigated 
were a "half station" approximation in which the finite difference approximations are 
made before expanding the derivatives of function products and a "whole station" approxi- 
mation in which derivatives of function products are expanded first; both approximations 
are in use. For the same number of stations, the average error in calculated deflection 
resulting from use of half station difference approximations was found to be always less 
than the error which would result from the use of whole station difference approxima- 
tions. The method was also applied to a square membrane subjected to a single Fourier 
type loading and a simple expression was obtained for the number of finite difference 
spaces required per Fourier wave length to keep finite difference results within a given 
percent error. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 24, 1967, 

124-08-06-29-23. 
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APPENDIX A 


BEAM WITH NONUNIFORM STIFFNESS 


As another example which illustrates the procedure described in this paper, con- 
sider a simply supported beam of unit length with nonuniform bending stiffness denoted 
by g, subjected to a uniformly distributed load of unit magnitude. The well-known dif- 
ferential equation governing the lateral deflection y of the beam is 

fe")" = 1 (Al) 

where primes denote differentiation with respect to x and variables are nondimension- 
alized to make the length of the beam, the bending stiffness at the left end, and the load 
each equal to one. The boundary conditions are 

y(x 0 ) = 0 y(x 0 + 1) = 0 

y”(x 0 )=° y"( x o + 1) = 0 

The left-hand side of equation (Al) is approximated by either the half station or 
whole station finite difference approximations for stations of equal spacing h. There- 
fore, from the half station approximation 

(gy'V - 1 = ^{fey")i-i ■ 2 fey”)i + fey'Vi] - 1 

= ^[ g i-l y i-2 " 2 ( g i-l + g i) y i-l + ( g i-l + 4g i + g i+l) y i 

" 2 ( g i + g i + l) y i + l + g i+l y i+2] - 1 < A2 > 

and from the whole station approximation 
(gy IV + 2g’y'" + g "y'Oj - 1 = 7 j(y i _2 * 4y i _ 1 + 6 yi - 4y i+1 + y. +2 ) 

+ + 2y *“l " ^+1 + y i+2) + ^"( y i-l " 2y i + y i + l) ' 1 


1 _ 

h4 


- h gi') y i-2 + (" 4g i + 2hg i' + h2g i") y i-l + ( 6g i ■ 2h2 gi")yi 
( _4g i - 2hg i’ + h2g i”)yi+i + ( g i + hg i')y i+ 2] - 1 ( A3 ) 
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APPENDIX A 


As before, expanding yj and gj about the ith point leads to the differential 
equation (6) where now 

L o(yi) = [g(x) i y i ’y'| 


B o (y±) = Vi 
B o 2 (yi) = yi" 

and for the half station approximation 




T , , gm 71 6i'yi v 7 „ IV em"' gi IVy i" 

L l( y i) “ + - 5 — + v> *i V* + —3— + — i T 

b iVi)*° 

iVi)‘f5 y i 




l ( y 0 6 ' 2 ' 12 D1 •' 1 3 12 

1, 


B x =.L 77 IV 


(A4) 


(A 5) 


J 


g^i 7111 gi'yi Vn 31 v! gi M, yi V 7 n , 

L 2 ( y i) 80 + 20 + 360 g i Yi + 12 + 144 g i y i 


gj v yi"’ gi^yi' 


60 360 


B n fVA = 0 


B 


} 2 ( y i) 

2 2 ( y i) 


— — — y . V! 
360 ^ 




(A6) 


and for the whole station approximation 


'i( y o = 


g^ gi' y i v ei" y i IV 


12 


Bj (yi) = 0 


B, 2 /y \ = — y.IV 

1 12 y ! 


J 


(A7) 
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APPENDIX A 


e . v .VIII e .- v .VII e .,, vi 

l 2 ( y i) = — + — + 


80 


20 


360 


*2 ( y 0 


B 2 2 (^i) 


.'l = 0 


J_ V.VI 
360 1 


J 


(A8) 


If solutions to equations (6) and (7), together with equations (A4) and either equa- 
tions (A5), (A6), . . ., or equations (A7), (A8), . . ., are again taken in the form of equa- 
tion (8), the series of simpler equations (11), (12), and (13) are again obtained (with 

Pi = 1) . Since the beam equation is fourth order rather than second, the boundary con- 

1 9 2 

dition of zero bending moment leads to specification of and B 2 . 

Results have been obtained for 


g(x) = X 11 


(n = 2, 3, 4) 
(11x12) 


for both the half station and whole station approximations of the derivatives. The error 
terms for both deflections and curvatures are shown in figure 2 for g(x) = x^ corre- 
sponding to a linearly tapered beam. An overall measure of the relative error in the 
half and whole station approximations is given in figure 3 for all three examples. The 
analytical error results for both deflection and curvature also agree with numerical 
error calculations within 1 percent for h less than about 0.1. 
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APPENDIX B 


CIRCULAR PLATE 


Clamped Circular Plate Under Uniform Load 

As another example of a second- order equation, consider the axisymmetric bending 
behavior of a clamped circular plate of radius 1.0 subjected to a uniformly distributed 
load of magnitude 2. If the plate has constant thickness and appropriate nondimensional 
variables are used, its behavior is governed by a second- order differential equation of 
the form , 

[!(*»] = -x (Bl) 

where x is the radial distance from the center and where 0 represents the slope of 
the plate. For a clamped plate the boundary conditions are 0 = 0 at x = 0 and x = 1. 

The two finite difference patterns for equation (Bl) are as follows: 

For the half station approximation, 


i(x0)’l +x i = — c 

L x Ji 


x i-l , 

x 1 ^i-1 l x 
1_ 2 




i_i “i+i 

1 2 1 + 2 , 


4 


+ Xj = 0 


and for the whole station approximation, 


4>i ^i 

H , 1 X 


0i + x: 9 + x i = 

x i x^ h 2 




+ x, = 0 


The differential operators Lg, Lj, and L2 in equation (6) are given by 

B o(4>i) = <t> i 


For the half station approximation, 

Ll (#i) = 


*i IV *i’" 4>l" $1 <f>i 

4- 




12 6x i 4xi 2 4xi 3 4xj 4 




! iW = 0 


J 


(B2) 


(B3) 


(B4) 


(B5) 
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T > <h n , *i v <h IV . *r *i" , *i 

^2 {jl) + i r>r> o + 7T 7 + 


*i 


360 120xi 48x^ 24x^ 16xj^ 16x^5 16xi® 


B 2 (^i) = 0 

and for the whole station approximation, 


■'n 


L, (* t ) = 


+ ■h'" 


12 6xj 


Bi(0i) = O 

v ^ <h v 

L 2 W = leo" + 120 Xi 

B 2 (0i) = O 




(B6) 


(B7) 


(B8) 


Results for the average principal error terms in the slope and in the numerically obtained 
second derivative obtained by using the previously described technique are presented in 
figure 3. 


Simply Supported Annular Plate Under Uniform Load 

The axisymmetric bending behavior of a circular plate is also governed by the fol- 
lowing fourth- order equation: 



where y represents the deflection of the plate. Results are obtained for a clamped 
plate annulus having an internal radius of 1 and an external radius of 2 and subjected to 
a uniformly distributed load of magnitude 1. The boundary conditions are y = 0 and 
y" = 0 at x = 1 and x = 2, respectively. 

The finite difference approximations to equation (B9) follow, and the results for the 
average finite difference error are given in figure 3. For the half station approximation, 
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The differential operators Lq, Lj, and L 2 in equation (6) are given by 


L 0(yi) [-(xy’)' 


B 0 ( y i) = y i 

B 0 2 ( y i) =y i' 

For the half station approximation, 


x _ x i y i IV . y i V y i IV 2 y i'" y i" y i’ 
L i( y 0 ~—r~ + ~9, — ~ + — ~ + 


3x i 3 Xi 2 Xi 3 Xj 


B 1 ( y i) = ° 


B i 2 ( y i) = ii y i IV 


j 


B 


^2 ( y i) : 
2 1 ( y i) 


x^i^ 11 y^ 1 


80 


2y i vi 2y. V y™ 2y.’" y." y.’ 

20 45 Xi 15x^2 3 Xi 3 3 Xi 4 x t 5 Xj 6 


= 0 


B 2 (y^ = — i— y.VI 

2 V V 360 


and for the whole station approximation, 


y i ¥i V V y i' 


'i (n) 


,, = 0 


B , 2 (m =I v ,IV 
1 (, y i) 12 y i 


J 


(B12) 


(B13) 


(B14) 


(B15) 
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yVm 13y.vn yV* yi v 

L 2( y i) 80 + 960xi 360xi 2 + 1920x i 3 

B 2 1 (y i ) = ° 

rx 2 , v \ = 1 y VI 
B 2 ( y i) 360 y i 


(B16) 


20 



APPENDIX C 


DEFLECTIONS OF A MEMBRANE 


As an example of the application of the method to a partial differential equation, 
consider a square membrane subjected to a unit sinusoidal loading and supported on all 
edges. The differential equation governing the membrane can be written as 


V 2 y = sin mirx sin nirz 

y = o 


(on boundary) 


(Cl) 


where the length of the side is 1 and where y is an appropriate nondimensional deflec- 
tion. Let equation (Cl) be approximated by finite differences with equal mesh spacing h 
in both directions. Two finite difference patterns which are often used to approximate 
V y are considered in this example. These operators are presented in symbolic form 
with their Taylor series expansions as follows: 

Standard approximation 


1 


1 

-4 

1 



(C2) 


"Refined” approximation 


1 

6h 2 


1 

4 


1 


4 

-20 

4 


1 

4 

1 


y = V 2 y + h 2 



4 - . . . 


(C3) 


where the subscripts denote partial differentiation with respect to the indicated variables. 
The refined approximation is denoted as such because it utilizes more node points than 
the standard approximation and for the special case for which the loading is linear (i.e., 
V^y es 0) has an order of error of h4 


The differential operators in equation (6) become 

L o(yi) = vl yii 

B o(yi) = yi 


(C4) 
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and for the standard approximation (C2) 

L l( y i) - y xxxx + y zzzz 

B l( y i) = 0 

and for the refined approximation (C3) 

L 1 (Yi) = v<r Yi 

B l( y i) = 0 

The finite difference solutions for the deflection at the ith point obtained by the 
perturbation method give 


(C5) 


(C6) 


Y 0 = 


sin m?rx sin n?rz 


^2^2 + n 2J 

and the principal error terms for the standard approximation 

,4\ 


(C7) 


1 + 




Yl = 


nr 


sin niTrx sin nirz 


(C8) 


12 1 + 


m‘ 


and the refined approximation 


Yi 


12 


sin mirx sin mrz 


(C 9) 


Although the refined approximation given by equation (C3) might be considered to be the 
better approximation, the error term shows that it is, in fact, less accurate for this 
problem. This result holds for all finite values of m and n; however, for m » n 
the errors in the two methods become essentially the same. 

Sample calculations were carried out to obtain actual numerical solutions and to 
compare them with the exact solution as well as with yj obtained by the perturbation 
method. Results were obtained for several values of h and m and n for both 
approximations and substantiate the greater accuracy of the standard approximation for 
this problem. The results are not shown but sample calculations for the dimensionless 
center deflection with h = 1/4 and m = n = 1 give 0.0533 by the standard approxima- 
tion and 0.0561 by the refined approximation; the exact answer is 0.0506. With h = 1/4 
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agreement was also obtained between the numerical results and the perturbation answers 
to three digits. This agreement could be improved to 5 digits if the h^ order error 
term Y 2 was included. 

Some practical assessment of the required number of stations to give a certain 
percentage error is also possible if equation (8) is written for the deflection as 

y i = Y 0 ^l + h 2 |i+ . . j (CIO) 

Let the principal error e be denoted 

e = h 2 (Cil) 

Y 0 

so that, for example a maximum error of 10 percent requires that h be chosen such 
that e < 0.1. For the standard approximation the principal error is 


e 




(C12) 


Since 1/ m is the length of a displacement Fourier wave, — ^ = N is the number of 
finite difference increments per wave length. Equation (C12) gives 


N = 


1 

mh 



(C13) 


which is fairly insensitive to n/m if m = n and which becomes for either m » n 
or n = m 


N~-!L= (C14) 

\fl2e 

For example, for a finite difference error of not more than 10 percent, N ~ 2.9. This 
means that approximately three finite difference spaces are required within the smallest 
Fourier wave length in order to obtain a 10-percent accuracy. For a 1-percent accuracy, 
9.1 spaces are required. Similar developments for the refined approximation give 


N « 



(Cl 5) 
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For m = n, equation (Cl 5) gives 



(Cl 6) 


or N « 4.05 for a 10-percent error and N « 12.8 for a 1-percent error. For m » n, 
the error is the same as that for the standard approximation. 
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x 



(a) f(x) = l/x3. 

Figure 1.- Finite difference error in deflection and curvature for a uniformly loaded string with nonuniform tension fix). 
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(b) f(x> = 1 + x3. 
Figure L- Concluded. 
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Figure 2.- Finite difference error in deflection and curvature of a uniformly loaded simply supported beam with nonuniform stiffness 

g(x) = > 
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(b) Second derivative. 

Figure 3.- Ratio of root-mean-square errors for half and whole station approximations. 
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